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velocity. The result is not consistent with the recent quasi-static theory of inelastic 
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1 Introduction 



The collision of particles with the internal degrees of freedom is inelastic in 
general. The inelastic collisions are abundant in nature (Goldsmith, 1960). Ex- 
amples can be seen in collisions of atoms, molecules, elastic materials, balls in 
sports, and so on. The study of inelastic collisions will be able to be widely 
accepted as one of fundamental subjects in physics, because they are almost 
always discussed in textbooks of elementary classical mechanics. 

Physicists realize that inelastic collisions can be a fashionable subject in physics 
from recent extensive interest in granular materials(Kadanoff, 1999; de Gennes, 
1999). In fact, granules consist of macroscopic dissipative particles. Therefore, 
the decision of interaction among particles is obviously important. We be- 
lieve that static interactions among granular particles can be described by 
the theory of elasticity (Love, 1927; Landau et al., 1960; Johnson, 1985; Hills 
et al., 1993). For example, the normal compression may be described by the 
Hertzian contact force(Hertz, 1882) and the shear force may be represented by 
the Mindlin force (Mindlin, 1949). The dynamical part related to the dissipa- 
tion, however, cannot be described by any reliable physical theory. Thus, the 
distinct element method (Cundall & Struck, 1979) which is one of the most 
popular models to simulate collections of granular particles contains some 
dynamical undetermined parameters. In other words, to determine such the 
parameters is important for both granular physics and fundamental physics. 

The normal impact of macroscopic materials is characterized by the coefficient 
of restitution (COR) defined by 

e = -VrjVi, (1) 
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where Vi and Vr are the relative velocities of incoming and outgoing parti- 
cles respectively. COR e had been believed to be a material constant, since 
the classical experiment by Newton (1962). In general, however, experiments 
show that COR for three dimensional materials is not a constant even in ap- 
proximate sense but depends strongly on the impact velocity(Goldsmith, 1960; 
Sondergaard et al., 1990; Bridges et al., 1984; Supulver et al., 1995; Giese et 
al., 1996; Aspelmeier et al., 1998; Basile et al, 2000; Labous et al, 1995). 

The origin of the dissipation in inelastic coUisions is the transfer of the kinetic 
energy of the center of mass into the internal degrees of freedom during the 
impacts. Systematic theoretical investigations of the impact have begun with 
the paper by Kuwabara and Kono (1987). Taking into account the viscous 
motion among the internal degrees of freedom, they derived the equation of 
the macroscopic deformation. Later, BriUiantov et al. (1996) and Morgado 
and Oppenheim (1997) derived the identical equation to eq.(2). In particular, 
the derivation by Morgado and Oppenheim (1997) is based on the standard 
technique of nonequilibrium statistical mechanics to extract the slow mode 
among the fast many modes which can be regarded as the thermal reservoir 
with constant temperature (see Appendix). Furthermore, BriUiantov et al. 
(1996) compared their theoretical results with experimental results. Thus, the 
quasi-static theory has been accepted as reasonable one. 

On the other hand, Gerl and Zippelius (1999) performed the microscopic sim- 
ulation of the two-dimensional collision of an isothermal elastic disk with a 
wall. Their simulation is mainly based on the mode expansion of an elastic 
disk under the force free boundary condition. The distinct characteristic of 
their model is that they do not introduce any dissipative mechanism of their 
microscopic equation of motion. Then, they solve Hamilton's equation deter- 
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mined by the elastic field and the repulsive potential to represent the collision 
of two disks. Their results show that COR decreases with the impact velocity, 
which strongly depends on Poisson's ratio. For high velocity of the impact they 
demonstrate the macroscopic deformation left after the collision is over. The 
relation between the quasi-static theory of impact (Kuwabara & Kono, 1987; 
Brilliantov et al., 1996; Morgado & Oppenheim, 1997) and their microscopic 
simulation (Gerl & Zippehus , 1999) is not trivial, because the energy transfer 
during the impact is not explicitly included in the quasi-static theory. Thus, 
we have to clarify the relation between two typical approaches. 

In this paper, we will perform the microscopic simulation of the impact of a two 
dimensional elastic disk with a wall. We introduce two methods of simulation: 
one is based on the lattice model (model A) and another is a continuum 
model (model B) which is identical to that by Gerl and Zippelius (1999). 
Both models do not include any dissipation explicitly. Thus, we regard inelastic 
collisions take place only from the transfer of modes of oscillation. Through 
our simulation, we will demonstrate that the elastic models do not recover 
the results predicted by the quasi-static theories in the low impact velocity 
(Kuwabara & Kono 1987; BriUiantov et al., 1996; Morgado & Oppenheim, 
1997; Schwager et al., 1998; Ramirez et al, 1999). 

The organization of this paper is as follows. In the next section, we will briefly 
review the outline of quasi-static theory (Kuwabara & Kono, 1987; Brilliantov 
et al., 1996; Morgado et al., 1997). In section 3, we will explain model A and 
model B which is equivalent to the model by Gerl and Zippelius (1999) of our 
simulation. In section 4, we will show the result of our simulation and discuss 
the validity of quasi-static theory. In section 5, we will discuss our results. 
In section 6, we will summarize our result. In Appendix, we summarize the 
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outline of the quasi-static theory by Morgado and Oppenheim (1997). Note 
that parts of this paper has been pubhshed as separated papers (Hayakawa 
& Kuninaka., 2001a; Hayakawa & Kuninaka, 2001b; Kuninaka & Hayakawa, 
2001). 



2 QuEisi-Static Theory: Review 



In this section, we briefly explain the outline of the quasi-static theory. One 
purpose of this section is to summarize the two-dimensional version of quasi- 
static theory which may not be mentioned in any articles explicitly. 

At first, let us summarize the three dimensional result, in which the equation 
of the macroscopic deformation is given by 

h = -kh{h^'^ + AhVhh) (2) 

in a coUision of two spheres. For the collision of two identical spheres the 
macroscopic deformation h is given by /i = 2R — |ri — with the radius 
R and the position of the center of the mass Fj of i th particle, h and h are 
respectively dh/dt and d'^h/dt^. kh in eq.(2) is written as 

where Mq, Yq and a are the mass of the sphere, (three dimensional) Young's 
modulus, and Poisson's ratio, respectively. In eq.(2) Ah is a constant, which 
may be a function of viscous parameters (Brilliantov et al., 1996). The first 
term of the right hand side in eq.(2) represents the Hertzian contact force 
(Love, 1927; Landau et al., 1960; Johnson, 1985; Hertz, 1882) and the second 
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term is the dissipation due to the internal motion. 

The simplest derivation of eq.(2) is that by BriUiantov et al. (1986), though 
we also check its validity by the alternative methods. Taking into account 
the limitation of the length of this paper, we follow the argument by them. 
The outline of the derivation by Morgado and Oppenheim (1997) is shown in 
Appendix. 

The static stress tensor in a two-dimensional linear elastic material can be 
represented by 



where /j, and K are respectively the shear modulus and the bulk modulus, and 
Uij is given by 



with the displacement field Ui. 

The two dimensional Hertzian contact law (Johnson, 1985; Gerl et al., 1999) 
is given by the relation between the macroscopic deformation of the center of 
mass h and the elastic force F^i as 



where Y is (two-dimensional) Young's modulus. Note that F^i and Y do not 
have dimension of the force and Young's modulus, because these are two di- 
mensional variables which are the ones per unit length along the third axis. 

Equation (6) can be derived from the stress tensor (4) with the standard treat- 



a'^'^'-'ij = 2ii{uij - SijUii/2) + KSijUii 



(4) 




(5) 




(6) 
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ment of linear elastic theory. Note that h satisfies h — R — yo with the position 
of the center of mass yo (Gerl & Zippehus, 1999). 

For small dissipation, as in Landau & Lifshitz (1960), the dissipative stress 
tensor due to the viscous motion among internal motions is given by 

a^^^'^j = 2r)i{uij - 5ijUii/2) + r)25ijuii, (7) 

where iiij is the time derivative of u^, rji {i — 1,2) is the viscous constant. 

Brilliantov el al. (1996) assumed that the velocity of deformation field is gov- 
erned by the macroscopic deformation, i.e., iii ~ h{dui/dh). Since in the limit 
of — > we may replace eq.(6) by F^i ~ —TrYh/ In^iR/h) (Gerl & Zippehus, 
1999). Thus, with the aid of the assumption by BriUiantov et al. (1996), (4) 
and (7), it is easy to derive the two dimensional version of quasi-static theory 



as 

nYh nYh 

where A is not an important constant. This result can be derived by various 
other methods. In section 4, we will compare the result of our simulation with 
eq.(8). 



3 Our Models 



Let us explain the details of our models to simulate collisions between two 
identical disks whose radius it! by the method of the mirror image. In both 
models, the wall exists at y = 0, and the center of mass of the disk keeps 
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the position at x — 0. The disk approaches from the region y > before 
rebounding from the wall. 

3.1 Model A 

The disk in model A consists of some mass points (with the mass m) on a 
triangular lattice. All the mass points are connected with linear springs with 
spring constant k. In the limit of a large number of mass points, this disk 
corresponds to the continuum circular disk with Young's modulus Y — 2k/\/3 
and Poisson's ratio 1/3 (Hoover, 1991). The position of each mass point of 
model A is governed by the following equation: 

= -^T.{do - \r, - r,|)^^^^ + eyaoVoe-^y^ (9) 

i=l nP ^il 

where do is the lattice constant, is the position of the nearest neighbor mass 
points of Tp, m is the mass of the mass points, is the y coordinate of r^, and 
By is the unit vector in the y direction. Note that the directional projection of 
the linear spring force in model A can cause the nonlinear deformation. The 
wall potential is given by Voe~"°^ , where Vq = mc'^aodo/2 with c = \Jy/ p 
and the density p. We adopt Oq = 100/do for the most of simulations, but 
we also adopt the result of ao = 25/do — 500/i? to obtain Fig.2, though the 
result is almost identical to that for a = lOO/do- The exponential interaction 
between the disk and the wall is introduced to simulate a collision between 
two identical disks. Actually, in the limit of oq — > oo, the exponential potential 
can be regarded as a potential of the mirror image. Thus, for later calculation, 
we analyze the case for large aodo. The number of mass points is fixed at 1459 
in model A, since the rough evaluation of convergence of the results has been 
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checked in this model. 



3.2 Model B 



In this subsection, we introduce model B which is originally proposed by Gerl 
and Zippelius (1999). Although the details of this model can be found in their 
paper, we present a short description of this model to understand the setup 
of our simulation. 

Gerl and Zippelius (1999) analyze Hamilton's equation to simulate collisions 
of a disk with the radius R as; 

dH OH 

under the Hamiltonian 

^ = 2M + + 2^<^^y + ^1 / #e--^(^'*). (11) 



n,l 



-7r/2 



Here M is the (two-dimensional) mass of an elastic disk, and Qn,i is the ex- 
pansion coefficient of the 2D elastic deformation field in the polar coordinate 

{ur{r,(t)),u^{r,4))) = 5^Q„,,«''(r)cosn(/),ii;''(r)sinn0), (12) 



1 nU \-n A dJn{kn^ir) Jn{k' , Jn{k' 

where u]!:'\r)R = An,i +nBn,i ■ — and m^' {r)R = -nA^r 



dr r r 

^ ^^^M^^n^ill -y^^ith the radius of the disk and the Bessel function of the n— th 
dr 

order Jn{x). Here k'^ i — kn,iyj2{l + a)/{l — cr^) and kn,i is the solution of 

(1 - a')(l - n^)KK'^Jn-i{K)Jn-i{K') + k^[k^ - 2n{n + 1)(1 - a)]Jn{K)Jn{K') 
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+(1 - a)[K^ - (1 - (7)(1 - n^)n][KJn-iiK)JniK') + J„_i ^(«)] = (13) 

with Poisson's ratio a, k — kn^iR and k' — k'^^iR, which is given by the force 
free boundary condition of tlic disk: 

(T,^(i?,0) = O (14) 



Thus, for fixed n there are infinitely many solutions kn,i and ujn,i = kn,iyY/ {p(l — 
numbered by / = 0, 1, • • • , oo. An^i and B^^i are determined by 

,AdJjk',R) JJk',R)^ , , 

+nS„(l - a)[- "^^'^ ^ - 1 = (15) 

and /o^ drr{it"''^ + it^''^} = -R^. -Pn,z is the canonical momentum. y(0, t) is the 
shape of the elastic disk in polar coordinates; 

Z/(0, ^) = Z/o(^) + X! Qn,i{Cn,i cos(n0) COS - ^n,/ sin(n0) sin 0) (16) 



with the position of the center of mass yo(t) and constants Cn,i and Sn,i deter- 
mined by the maximal radial and tangential displacement at the edge of the 
disk as Cn,i — u^'\R) and Sn^i — u^'\R). M is the mass of the disk, and the 
momentum of the center of the mass Pq = Myo satisfies Pq = —{dH/dyo) , Vq 
and a are parameters to express the strength of the wall potential. 

For the simulation of a pair of identical disks, they have confirmed that the 
result with finite can be extrapolated to the result of ao — > oo by taking 
into account finite ao effect in proportion to l/(aoi?). Similarly, the result 
with finite number of modes N should be extrapolated with the correction in 

proportion to 1/\/N. Since they have already checked such the tendencies, we 
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only adopt N = 1189 {n < 50 and Kn < 50)or N = 437 {n < 30 and k„ < 30), 
Vi = Mc^aoR/2 and ao = 500/ R. 



3.3 Parameters in both models 

For the comparison between two different models, we only simulate the case of 
Poisson's ratio a — 1/3. The numerical integration scheme for model A is the 
classical fourth order Runge-Kutta method with At = 1.6 x 10~^^m/ n. Parts 
of the calculation in model A has been checked by the fourth order symplectic 
integral method with = 5.0 x 10~^^m/ k., and no differences in results of 
two methods can be found. For model B, we adopt the fourth order symplectic 
integral method with At — 5.0 x 10~^i?/c. In both models, we have checked 
for conservation of the total energy. 

We also investigate the impact with finite temperature. The temperature is 
introduced as follows: In model A, we prepare the Maxwellian for the initial 
velocity distribution of mass points, where the positions of all mass points are 
located at their equilibrium positions. From the variance of the Maxwellian we 
can introduce the temperature as a parameter. To perform the simulation, we 
prepare 10 independent samples obeying Maxwelhan with the aid of normal 
random number. In model B, we prepare samples which satisfies Gibbs states. 
Namely, \/MuJn,iQn,i/ and Pn,i/y/'2.M obey the normal random number 
with the variance (temperature) T. In model B, we prepare many samples 
(120 or 20) to simulate systems at finite T. 

The summary of differences between model A and B is as follows: (i) All of 
the mass points in model A interact with the wall but, in model B, only the 
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exterior boundary has the influence of the potential as in eq.(ll). We have 
replaced the original model A by a model in which only mass points on the 
boundary can interact with the wall, but we cannot find significant differences 
in the results of our simulation in both discrete models, (ii) Model A can 
have nonlinear deformations, but model B is based on the theory of linear 
elasticity, (iii) Model A can express some plastic deformations, but model B 
cannot. This effect will be discussed in section 6. (iv) Model A has six fold 
symmetry whereas model B has only rotational symmetry, (v) The force free 
boundary condition (14) is assumed in model B but may not be appropriate 
for actual situations. Model A does not include such the condition. 

4 Results 

Now, let us explain the details of the result of our simulation. In the first 
subsection, we will introduce the result at T = and in the second subsection, 
we will show the result at finite T. 

4-1 Simulation atT — 

At first, we carry out the simulation of model A and model B with the initial 
condition at T = (ie. no internal motion). Figure 2 is the plot of the COR 
against the impact velocity for both model A and model B. For model A, we 
have adopted the fourth order Runge-Kutta method. To ehminate the effect 
of six fold symmetry of model A, we average 12 data as a function 9 of the 
initial orientations of the disk i.e. 9 — 7rn/72 with n — 1, 2, 3, • • • , 12 with 
ao = 25/ do = 500/ R for = 1459. We also investigate the case that only 
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mass points at the boundary can interact with the wall for small Vi but their 
results do not have any visible difference from the original model A. It is 
obvious that there is no plastic deformations for Vi < 0.2c. 

For model B, we show the results of 437 modes and 1189 modes which clearly 
demonstrates the convergence of the result for the number of modes. When 
impact velocity Vi is larger than 0.1c with c = ^Y/ p, the value of COR of 
model A is almost identical to that of model B. Each line decreases smoothly 
as impact velocity increases. 

At present, we do not know the reason why the significant difference between 
the two models exists at low impact velocity. It is difficult to imagine that 
occurrence of nonlinear deformations during the impact of model A causes 
the difference because the deformation is smaller when Vi is smaller. 

Second, we investigate the force acting on the center of mass of the disk 
caused by the interaction with the wall in model B. In the limit of ^ we 
expect that the Hertzian contact theory can be used(Landau & Lifshitz, 1960; 
Johnson, 1985; Gerl & Zippelius, 1999). The small amount of transfer from 
the translational motion to the internal motion is the macroscopic dissipation. 
Thus, we can check whether the quasi-static approaches (Kuwabara & Kono, 
1987; Brilliantov et al., 1996; Morgado & Oppenheim, 1997) or our elastic 
simulation can be used in slow impact situations. 

If h is given, we can calculate the elastic force by solving eq.(6) numerically. 
Figure 3 is the comparison with our simulation in model B (1189 modes) and 
the Hertzian contact theory (6) which is given by the solid lines. The result 
of our simulation at the impact velocity Vi — 0.01c shows the hysteresis as 
suggested in the simulation at Vi = 0.1c (Gerl & Zippelius, 1999). This means 
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the compression and rebound are not symmetric. The hysteresis curve is still 
self-similar even at Vi — 0.04c but the loop becomes noisy at Vi — 0.1c. 

For very low impact velocity Vi — 0.001c, the hysteresis loop almost disappears 
and the total force observed in our simulation is almost a linear function of h 
which deviates from the one predicted by both the Hertzian contact theory and 
the quasi-static theory (8). In particular, the turning point which corresponds 
to the point of the largest in Fig. 3(b) is apart from the Hertzian curve 
(the sohd hue). This deviation is in clear contrast to the quasi static theory, 
because the dissipative force in the theory in eqs.(2) and (8) must be zero at 
the turning point which h — should satisfy. This tendency is invariant even 
for the simulation of model A, though the data becomes noisy. The hnearity 
of the total repulsion force is not surprising, because e"""^*^"^'*^ in the potential 
term in eq.(ll) can be expanded in a series of Qn,i for very slow impact. 

The result may suggest that our elastic models do not recover the Hertzian 
contact theory in the quasi-static limit. To check the tendency, we investigate 
whether any static state can be achieved in our models in the compression. 
Figure 4 is the time evolution of the center of mass in the simulation of model 
B, where the strength of dimensionless external field is g = O.Olc^/R. We 
observe that an undamped harmonic oscillation of the center of mass in the 
simulation after the first deformation. This oscillation is stable because the 
energy of oscillation is not enough to overcome finite energy gap between en- 
ergy levels. Thus, the center of mass keeps the oscillation as the motion in 
the ground state. We note that Fig.4 is the result of the simulation at finite 
temperature in which the mode transfer is enhanced. Nevertheless, the cen- 
ter of mass keeps the harmonic oscillation. This tendency can be observed 
in model A, too. Even when we introduce the randomness in the coupling 
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in model A, the oscillation is undamped. Thus, both of elastic models can- 
not reach any cquiUbrium steady state as is assumed in the Hertzian contact 
theory. This result indicates that the elastic models are not appropriate to 
describe quasi-static situations for i>j/c <^ 1. Note that the introduction of 
nonlinear deformation may not be enough, because as we can sec in Fig. 3 (b) 
the deformation is very small for slow impact. Thus, it is difficult to imagine 
the impact produces nonlinear deformations. To reach an equilibrium state, 
thus, we need to introduce some microscopic dissipative mechanism. 

However, the validity of the contact time r in the impact evaluated as r ~ 
{t:R/ c)^\n{Ac/ Vi) by the quasi-static theory (Gerl & Zippelius, 1999) has been 
confirmed by the results of our simulation of model A (Fig. 5). Thus, our elastic 
model can be valid in the impact with the intermediate speed. 



4-2 Simulation at finite T 

Now, let us show the results of our simulation at finite T. The thermal velocity 
Vth — \Jt/M causes significant differences from those at T = in both low 
and large impact velocities. In this sense, we have much room to study this 
process at finite T systematically. 

For small impact velocity, i.e. if the effect oivth is not negligible, the fiuctuation 
of COR at finite T becomes large, while the average is almost independent 
of temperature as in Figs. 6 and 7, where the results are obtained from the 
average of 120 independent samples. In some trials at high temperature, thus, 
COR becomes larger than 1, though the average is less than 1. Of course, 
for such the high temperature, it is impossible to control the actual speed of 
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impact. 



For large impact velocity, Vi ^ Vth, we do not observe any definite temperature 
effect in model B but wc find drastic decrease of COR in model A. It seems 
that COR can be on a universal curve when the impact velocity is scaled by the 
critical velocity above which the COR decreases abruptly (Fig.8). The relation 
between the critical velocity and the initial temperature at the intermediate 
impact velocities is shown in Fig. 9. The critical velocity seems to obey a linear 
function of T, though the data is not on the function for both slow and fast 
impacts. 

5 Alternative Quasi-Static Theory: The Effect of Temperature Gra- 
dient 

In this section, let us discuss new aspects of the quasi-static theory. As in 
section 2, the conventional quasi-static theories (Kuwabara & Kono, 1987; 
Brilliantov et al., 1996) consider the effect of internal friction. Similarly, the 
Langevin approach (Morgado & Oppenheim, 1997) gives the identical result 
to that by conventional one. In both approaches, it is assumed that the tem- 
perature in disks is uniform. However, this assumption is not accurate. It is 
known that the rise of temperature is proportional to the divergence of elas- 
tic deformation (Landau & Lifshitz, 1960). Thus, the temperature cannot be 
uniform. 

In this section, we will evaluate the dissipation rate due to the thermal diffu- 
sion and show that the contribution of this term is dominant in quasi-static 
situations. The result may not be complete but meaningful to indicate the 
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importance of the thermal diffusion. 



In a quasi-static coUision, the compression is proceeded in an adiabatic process. 
The adiabatic condition is written as So(T) + Kauu = So(Tq), where 5*0, K, 
a and Tq are respectively the entropy (divided by the Boltzmann constant), 
the bulk modulus, the thermal expansion rate and the temperature without 
any deformation (Landau, 1960). From the expansion of the entropy around 
To we obtain 

T -Tq^ Uii = T^(C; - Ct)Uii, (17) 

Op Op 

where Kad, Cp, Q and q are the bulk modulus in the adiabatic process, the heat 
capacity at constant pressure, the sound velocity of the longitudinal mode and 
the sound velocity of the tangential mode, respectively (Landau & Lifshitz, 
1960). To obtain the final expression we use the two-dimensional relations 
= Y/{2{1 - a)), ci = ,/Y/(Xr^ and q = ^jY/{2p{l + a)). There is 
the relation between the stress tensor and the divergence of deformation uu 



as 



1 — 0" 1 — O" 

Uii = y an = y [a^x + ayy) (18) 



in the two-dimensional elastic medium (Landau & Lifshitz, 1960). Substituting 
(18) into (17) we obtain 

T-To = ^au. (19) 

Thus, if (Tjj is a function of the position, the temperature field is not uniform, 
which is contrast to the assumption in previous quasi-static theory. 

It is known that the thermal diffusion causes energy dissipation. The dissipa- 
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tion rate is given by 

E--yJ d'v{VT)\ (20) 

where ht is the thermal conductivity (Landau & Lifshitz, 1960). The integra- 
tion in (20) is performed in all region of elastic disks. Thus, from (19) and 
(20), the energy dissipation which is not included in previous treatments is 
need to be considered. 

Now, let us evaluate the integral (20). For this purpose, we use the exact 
solution of two-dimensional Hertzian contact problem (Hills et al., 1993). The 
explicit stress tensor is given by 



S V i + X s 



V/TT^ S (1 + 52)3/2(^4 + ^2) 



"- = -^°IFTr)' 

where x — x/a and y — y/a are scaled by the contact radius a which is given 
by 



for the contact of two identical disks. Note that x and y are the position in 
the Cartesian coordinate whose origin is the center of the contact area (see 
Fig. 10). po in (21) is given by 

Po = (23) 
na 

where s in (21) is 

I {-(1 f) + ^(l-x2-y2)2 + 4^2| . (24) 



2 
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Prom (21) we obtain an 



Thus, the dissipation rate (20) can be calculated in principle. 

Note that the numerical integration of (20) is not easy, because (i) the explicit 
expression is too complicated, (ii) the boundary is modified by the compres- 
sion, and (iii) the parameter R = R/a is important and is a function of the 
impact velocity. Thus, here, we present a rough analytical evaluation of (20) to 
capture the characteristics of this problem. We note that an becomes simple 
in some special situations. For example, an at x — which is on the axis of 
symmetry is given by (Hills et al., 1993) 



ai^ = au{0,y)^2po 



(26) 



On the other hand, the integral representation of ctj, 



(^ii = / «?7 — 2' (27) 



can be approximated by 

.r'--%/<i£V^=-% (28) 

—a 

far from x — 0. Here we use /"^ d^^/d^ — = 7ra^/2. 

For the evaluation of (20), we distinguish the inner part |x| < a from the outer 
part \x\ > a. In the inner region, we may replace an by a*". Thus, (VT)^ in 



19 



the inner region may be approximated by 

(v^)L.^(l-^)^ (29) 



In the outer region we may replace an by because such the approximation 
can be used in the most of regions in the quasi-static situation {R^ 1). Thus, 
(VT)^ in the outer region may be approximated by 

(vrjL - ^i^r (30) 

Of course, these assumptions cannot be used in general. In particular, near 
the edge |a;| = a the contribution is expected to be large. However, we believe 
that the evaluation under the simplified assumption may be useful as the first 
step of the analysis. 

In the inner region, the integrand is independent of x and the integrated region 
may be approximated as a square domain —a < x < a and < y < 2R. Thus, 
Ein in (20) can be evaluated as follows: Prom 

2R I — ^ 

I dt(l - ^jJ=r = 2 + 4^(1 + tan-(2^) ^ (31) 

with i? ^ 1, we obtain 

For the outer region, Eout in (21) is 



The evaluation of the outer region is more complicated, because the domain 
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can not be approximated by a simple rectangular domain. For the evaluation, 
we neglect the deformation of shape of the compressed disk. Thus, the shape 
is approximated by a hemi-circle as in Fig. 10. It is convenient to introduce 
the polar coordinate (f, 9) to evaluate (33). For a given angle 9 between x axis 
and OQ in Fig. 10 r is between f„jj„ = OP /a and fmax = OQ/a. We also 
introduce 9min which is the cutoff angle for 9. Taking into account ^ ^ 1 we 
can evaluate 

r^„^~2^sin^: r^^^ ~ ^— . (34) 

COSb* 



Since these evaluations are approximate, we need to introduce the lower cutoff 

of 9min by the consistency condition fmax{9min) > rmin{9min) = 1- Thus, 9min — 

a/2R. The upper cutoff of 9 is 9max = cot-i(l/2^) ~ 7r/2 - 1/2^. Thus, the 
integral in (33) can be evaluated as 



'f'max (^) 



(0) 

= ^ / ^^ + 7; / d9cos^9c^- ^^77, (35 

8R\J sin^9 2 J 8 2i? 8' ^ ^ 

"min "min 

where we use 



I 



' d9 1,1,^2 

= ^+cotH^ ~2i?+^ 36 

sm^^ 2R ^2R 3R 



with cot9 ~l/9-9/3 in the limit of ^ ^ 0, and 

/ d9cos'9^--X. (37) 
i 4 4i? ^ ' 
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Substituting (35) into (33) we obtain 



E, 



'out — 



32C/a2 ■ 



(38) 



From (32) and (38), the total dissipation rate E = Ein + Eout is given by 



E = -7o 



KTToYa^Fei 
C/{l-a^)R 



(39) 



where 



7o = 



TT^ + 128 - 327r 
2567r2 



(40) 



The result suggests that the dissipation rate by the thermal diffusion is dom- 
inant in quasi-static situations, because the force F^i appears in (39) exists 
even in the limit of zero impact velocity, while internal frictions considered in 
conventional quasi-static theory disappears in the limit of zero impact velocity. 

This result, however, predicts a singular behavior of COR. In fact, the rough 
evaluation of the total energy loss Ei^ss by heat diffusion during the impact is 
proportional to the impact velocity Vi, while the definition of COR by Eioss 
is Eioss = Mv^{l — e^)/2. Thus, COR may be singular for very small impact 
velocity. We need to consider another mechanism to remove such the singu- 
larity. We also need such the analysis for three dimensional situations where 
the stress field becomes simpler than that for two-dimensional cases (Hills et 
al. 1993). 
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6 Discussion 



We investigate what happens in the disk above the critical velocity and find 
the existence of plastic deformation of the disk (Fig. 11(a)). Actually, there 
are no energy differences between two configurations in Fig. 11(b) which can 
occur after the strong compression during the impact but cannot be released 
after the impact is over. It is well known that plastic deformation causes the 
drop of the COR (Johnson, 1985). 

6.1 Application of the Conventional Theory of Plastic Deformation to 2D 
Impacts 

Following the description by Johnson (1985), let us explain the dimensional 
analysis of the two-dimensional plastic deformation. From two-dimensional 
Hertzian law (6) we evaluate h ~ a'^/R (Johnson, 1985). The work for the 
compression of the disk W is W ^ (l/2)Mv| ~ fo* dhF^ ~ ff daa^/R^, 
where M and Vi are the mass of the disk and the impact velocity, respectively. 
h* and a* are respectively the maximal compression and and the maximal 
contact length. Here we neglect the logarithmic correction and unimportant 
numerical factors. Introducing the mean contact pressure during dynamical 
loading which satisfies Pd ~ F^i/a, W can be evaluated by 1^ ~ F^ih* ~ 
Pd{a*)^/R. From W ~ Mvf we can express a* ~ (MyfR/pd)^^^. 

Let us assume that the impact exceeds the yield pressure for the plastic de- 
formation. In such the case, the deformation during rebound is frozen. Thus, 
the work in a rebound is W ~ F*h* where F* is the maximal force during 
the impact. From h* ~ F*/Y and F* ~ p^a* we evaluate W ~ {pda*)^/Y. 
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Substituting the expression of Oq into the expression for W and W we obtain 
the COR as 

vf W ^ Y{Mvfy/^' ^ ' 

— 1/3 

Thus, we expect the law e ~ f j in the coUision of a plastic deformed disk. 
The three dimensional version of evaluation which gives e ~ v^^^'^ agrees well 
with the experiment (Johnson, 1985). 



6.2 Realistic Systems 



The actual plastic deformation is more complicated than what we modeled in 
this paper. For example, in the actual contact area a central region of perfect 
contact is surrounded by an annulus of imperfect contact. In actual situations, 
it is not easy to obtain a pure normal collision, because the rotation of disks 
is difficult to be suppressed and the wall is not perfectly flat. Thus, a little 
deviation of the collision angle causes the tangential stress in collisions. In 
the existence of tangential stress, we need to consider the effect of imperfect 
contact or partial slip in the outer region to get finite force at the corner of 
contact area. 

We also note that the actual materials are not uniform. They contain a lot 
of microcracks, and amorphous structure locally. Such the imperfection of the 
materials causes the local achievement of the yield of plastic deformation. 
Thus, the plastic deformation also occurs localized in contrast to the macro- 
scopic deformation in Fig. 11. 

Our finding is, however, something new, because (i) the decrease of COR is 
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excited by the temperature and (ii) COR decreases more rapidly like e ~ 
v~^'^ than that for the conventional plastic deformation e ~ in (41). 

The mechanism how to occur the plastic deformation is not clear at present 
including the linear law in Fig. 9. 

For future refinement of our model to describe plastic deformation, we need to 
introduce (i) the initial cracks, (ii) local deformation of lattices at the initial 
condition, (iii) the yield of local plastic deformation or non-Hookian effects of 
springs, and (iv) porosity distribution at the initial condition except for the 
introduction of the heat diffusion effects as introduced in section 5. Of course, 
to compare the simulation with experiments, we have to simulate the model 
in three dimensional situations. 

7 Conclusion 

We have numerically studied the impact of a two dimensional clastic disk with 
the wall with the aid of model A and model B. The result can be summarized 
as (i) The coefficient of restitution (COR) decreases with the impact velocity, 
(ii) The result of our simulation is not consistent with the result of the two- 
dimensional quasi-static theory. For large impact velocity, there is hysteresis 
in the deformation of the center of mass. For small velocity, there remains the 
inelastic force even at /i = 0. (iii) The effect of heat diffusion may be important 
for the small impact velocity, (iv) There are drastic effects of temperature 
in both small and large impact velocity, (v) In particular, for large impact 
velocity of model A, we have found the abrupt drop of COR above the critical 
impact velocity by the plastic deformation. The critical velocity of the plastic 
deformation seems to obey a simple linear function of temperature. 
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We believe that this prehminary report is meaningful to recognize that physi- 
cists have poor understanding of such the fundamental process of elementary 
mechanics. We hope that this paper will invite a lot of interest in the impact 
from various view points. We, at least, have a plan to study three dimensional 

impacts to clarify the relation among the microscopic simulation, experiments 
and the quasi-static elastic theory. 



Appendix 



A Langevin Approach to the Quasi-Static Theory 

In this Appendix let us introduce the derivation of quasi-static theory by 
Morgado and Oppenheim (1997). The characteristics of their derivation is to 
introduce 'thermal deformation' explicitly to assist the elastic deformation. 
Although they do not mention what the thermal deformation is, it is the com- 
plex combination of inelastic scattering of phonons, electrons , sound radiation 
into the air and any other mechanism which cannot be regarded as the elastic 
deformation. In this Appendix we introduce a simplified version of the deriva- 
tion of Langevin equation instead of using the original argument (Morgado & 
Oppenheim, 1997). Note that the argument in this Appendix is restricted to 
one in three dimensional systems. 

Thus, the position dj of i-th. atom can be written as dj = -|- Uj -|- Pj, 
where Rj, Uj, and Pj are respectively the equilibrium position of the atom i, 
the elastic deformation and the 'thermal' deformation. Here the terminology 
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of 'thermal' deformation means that the deformation cannot be controlled or 
the origin has not been specified from macroscopic point of view. Such the 
'thermal' deformations may be regarded as a not important variables which 
can be treated as random variables. Let us introduce = dj — Rj. In the local 
rule in this Appendix, suffices i,j represent atoms, and the Greek suffices 
such as a, (5 are components. In addition, we adopt Einstein's rule for suffices 
where the duplicated suffices mean the summation as ajja = Ha^a^a- For 
small deformation, can be written as 

I, ~ R, ■ — u(R,) + p,. (A.l) 

The potential among atoms can be approximated by a harmonic one near 
its equilibrium position. Thus, the harmonic potential for isotropic systems is 
given by 

t/o = fEe (A.2) 

Note that the original paper(Morgado & Oppenheim, 1997) does not assume 
isotropic form of Uq. Then Uq is replaced by a second order tensor. Substituting 
(A.l) into (A.2), eq.(A.2) becomes the combination of three terms: 

Uo = Uei + U^ + Uh (A.3) 

The first term of the right hand side of (A.3) is the elastic energy for the 
deformation which can be written as 

Uel = X){^'"aa(Ri)'"/3^(R'i) + IJ'Uap{'Ri)Uap{'Ri)} , (A.4) 

where uRaRp = + |^a7^7/3 Note that in the original paper (Morgado & 
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Oppenheim, 1997), the coefficient of Uapu-ys becomes a fourtli-order tensor. 
also includes a constant which is represented by a second-order tensor. 



Here, A and jji are Lame's elastic coefficients. The second term of right hand 
side of (A. 3) is given by 



which expresses the coupling between the elastic deformation and the thermal 
deformation. The third term of (A. 3), Uh = ^ PiaPja, is the energy of the 
thermal deformation. The contribution of this term is in general smaller than 
other terms. 

The collision of two elastic bodies consists of materials 1 and 2. The energy is 
the simple summation of the contribution of two materials. Let the center of 
mass of i-th particle r^'^ (i = 1, 2). The interaction during coUision appears if 



Morgado and Oppenheim (1997) assume that the slow mode, the motion of 
the center of mass, can be written by the solution of the Langevin equation. 
In the Langevin equation, the elastic force is regarded as a systematic force, 
while the force —VU^ plays a role of the fluctuating force. As in the general 
framework of the fluctuation-dissipation theorem, the friction coefficient C(^i2) 
in the Langevin equation is determined by the time correlation function of the 
fluctuating force as 




ri2 = 



r(i) -r(2)| < 2R. 



oo 




(A.5) 
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Here VU^{t) is retarded one of VU^ by time r. Thus, we can write 



where V12 = and the upper suffix {k) represents the particles 1,2. Here, 
Pj can be regarded as the thermal fluctuation as 



K 



(A.7) 



Thus, we obtain 



C(ri2) = E 4'=)<(V,24'^)(Vi2<)). 

I 1 



fe=l 



(A.8) 



Let us recall that Kramer's equation for many-body systems 



dt 



+ 



i=l 
2 



M 



i=l 



P{Xt,t) (A.9) 



is equivalent to the Langevin equation, where (3 = 1/T, Xt = {pj(t), rj(t)} 
{i = 1,2), C12 = Ci^u), Pkj = Pj - Pfe, Vp^. = Vp. - Vp,. fj-fe is the unit 
vector from the center of j th. particle to the center of fc-th. particle. Here we 
introduce the average < B >t as < B >t= J dXtBP{Xt, t) for any variable 
B. Note that we have the relation 



dP 



(A.IO) 



Introducing the relative coordinate ri2 = r2 — ri approximating that the 
potential Uq is approximated by Ugi we obtain 
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< Pl2 >t = - / dXtP{Xt,t){Vr,Uel - Vr,Uel) - J dXtCuhh " Pl2P{Xt,t) 
= < F^i >t -J^ < Cl2fl2fl2 • Pl2 >t ■ (A.ll) 



Note the contribution of the hnear momentum disappears from the integral 
by parts. The elastic force < ¥ei >— —ri2U'^i{h) is nothing but the Hertzian 
contact force, and Uei = ^h^^"^. Thus, the time evolution of the macroscopic 
deformation h {h satisfies h = 2R — ru-) is given by 



EK=-h,h^n-mh. (A.12) 



Therefore, to determine the dissipation is reduced to determination of the 
friction constant C(^)- 

Unfortunately, it is impossible to obtain the exact form of C(^)) because ( 
is determined from the complicated relations between the elastic deforma- 
tion and the thermal deformation. However, from the consideration of power 
counting of h it is not difficult to deduce how (^{h) depends on h. In fact, 
it is easy to show the scaling Mz(x) Uz{y/a:>L) = aUz{'x) in the Hertzian 
contact theory. Similarly we have Ma^(x) —>■ Ma^(-\/ax) = s/aua/si^) 
dua/si'x.) / dh — > duap{\/ayi) / adh — l/\/a{duai3{y^/dh). Prom the comparison 
between the elastic energy (A. 4) and C(^) in (A. 8), it is easy to understand 
that the key point is how Uap and dua/s/dh are scaled by a. From the discus- 
sion here the elastic energy is scaled as a^/^, and thus ({h) is scaled as a^/^. 
Thus, we finally obtain 

ah) = ^k'h'/' (A.13) 
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and 



fh = -lhh'/'-h"h, (A.14) 



where k" = j3k'/M and k' cannot be determined from this argument. This 
result agrees with the result by the viscous stress tensor (Kuwabara & Kono 
, 1987; Brilliantov et al., 1996). 

Note that the derivation is quite different from the previous one assumed the 
existence of viscous tensor. Both of derivation assumed that the temperature 
field is uniform, but this assumption is not correct in general. As discussed in 
the text, the rise of temperature is directly related to the compression. Since 
the compression is not uniform, the rise of temperature is not uniform. 

Two-dimensional quasi-static theory in eq. (8) can be derived by a parallel 
argument introduced in this Appendix. 
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Notation^ 



In this paper, we analyze two dimensional systems. As a result many of quantities 
have different dimensions from those for usual three dimensional ones. We also note 
that we do not introduce Boltzmann's constant in the calculation. Thus, T has the 
same dimension with E, and the entropy becomes dimensionless. 
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oo coefficient of a wall potential, 1/m 

a contact radius, m 

A,Ah constants, s 

c one dimensional sound velocity, m/s 

Q sound velocity of the longitudinal mode, m/s 

Cp heat capacity at constant pressure, dimensionless 

Ct sound velocity of the tangential mode, m/s 

do lattice constant, m 

e coefficient of restitution, dimensionless 

By the unit vector of y direction. 

E energy (2D), J/m 

Eioss energy loss during an impact (2D), J/m 

E the dissipation rate (2D), J/m- s 

Ein the dissipation rate in the internal region (2D), J/m- s 

Eout the dissipation rate in the outer region (2D), J/m- s 

Fei 2D elastic force, N/m 

Ftot 2D total force, N/m 
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g an external field, dimensionless 

h macroscopic deformation, m 

H 2D Hamiltonian, J/m 

Jn{x) n-th order Bessel function 

kh the constant in eq.(3) 

K 2D bulk modulus, N/m 

Kad 2D bulk modulus in the adiabatic process, N/m 

m mass of a mass point (2D), kg/m 

M mass of an elastic disk (2D), kg/m 

Mo mass of a sphere, kg 

Pn,i canonical momentum (2D), kg/s 

Po the constant in eq.(23) 

Qn,i canonical coordinate, m 

R radius of the disk, m 

R dimensionless radius (R/a) 

Fj position of mass point 

^mao; maximum distance from origin, dimensionless 
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fjnin minimum distance from origin, dimensionless 

f distance from origin, dimensionless 

s defined in eq.(22) 

5*0 dimensionless entropy, 

t time, s 

T temperature (2D), J/m 

To temperature without deformation, J 

Ui displacement field, m 

Uij strain tensor, dimensionless 

Vq coefficient of wall potential, J/m 

Vi impact velocity, m/s 

Vr rebound velocity, m/s 

Vth thermal velocity, m/s 

X dimensionless x coordinate of the position (x/a) 

Y Young's modulus (2D), N/m 

Yq Young's modulus (3D), N/m^ 

Ho position of center of mass, m 
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yi<P,t) 
y 



y coordinate of r^, m 

the shape of an elasic disk, m 

dimensionless y coordinate of the position (y/a) 



Greek letters 

a 

Vi 
To 

K 

Kt 
P 

a 



thermal expansion rate (2D), m/N 
1/T 

viscous constant (2D), N- s/m 
the dimensionless constant in eq.(40) 
spring constant, N/m 
thermal condutivity (2D), 1/s 
shear modulus (2D), N/m 
density (2D), kg/m^ 
Poisson's ratio, dimensionless 
2-dimensional static stress tensor 
stress tensor in the inner region 
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CTj""* stress tensor in the outer region 
At time step for simulation, s 
T contact time, s 
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Figure 1: A schematic figure of a disk used in model A. 

Figure 2: Coefficient of restitution for normal collision of the Model A and 
Model B as a function of impact velocity, where c = I P with Young's 
modulus Y and the density p. 437 and 1189 modes are chosen for model B. 
The error bar in model A represents the standard deviation of the data as a 
function of the initial orientation Q. 

Figure 3: The comparison of the Hertzian force in eq.(6) with our simulation 
at Vi — 0.01c (a) and Vi — O.OOlc(b) at T = in model B. Ftot is the total 
force originated from the interaction with a wall. 

Figure 4: The time evolution of the center of mass of the elastic disk under 
the compression by = O.Olc^/i? (model B with N=437). Here the dimen- 
sionless time is measured by Rjc and the position of CM. (center of mass) 
is measured by the diameter of the disk (2R). Simulation is performed at the 
finite temperature T — 10~®Mc^ and is averaged over 20 independent samples 
which start from initial condition satisfying the Gibbs distribution. 

Figure 5: The plot of contact time versus the impact velocity (model A). R 
represents the radius of the disk, in which R = 40,60,70 and 80 correspond to 
the number of mass points 5815, 13057, 17761 and 23233, respectively. The 
dash-dotted line is fitting curve based on the quasi-static theory. Here the line 
represents the fit of tc/R as 3.21758y^ln(4c/t'i) ~ TT^ln^ic/vi). 

Figure 6: The average shift of COR at finite temperature T = 10~^Mc^ as 
a function of the impact velocity in model B with N — 437. The dotted hne 
indicates one at e(T) = e(0). 

Figure 7: The standard deviation of COR a = J< (e- < e >y > at T = 
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10 ^Mc^ as a function of the impact velocity Vi via model B with N — 437. 

Figure 8: The relation between the coefficient of restitution and the impact 
velocity rescaled by the critical velocity for each temperature. Curves are 
plotted in the log-log scale. The temperature is scaled by Tq — mc^ with the 
mass of the mass points m. Note that the error bars are plotted only in the 
case T/Tq = 0.03 but are the same order even at other T (model A). 

Figure 9: The plot of the initial temperature and the critical velocity causing 
the plastic deformation. Vcr/c — a{T/TQ) +b is the fitting curve line from the 
data between T/To = 0.02 and 0.05 (model A). 

Figure 10: The configuration of a compressed disk. The origin is O. The angle 
between OP and x axis is 9. The inner region is inside two vertical dashed 
line, while the outer region is outside the central region. The length of OS is 
equal to a. The radius of (undeformed) disk is R. 

Figure 11: (a) Plastic deformation of model A with Vi — 0.22c at T = 0.03mc^. 
The solid circle represents the initial circle. The points in a circle are positions 
of the mass points after the collision. The deformation is asymmetric because 
of the velocity distribution at the initial stage, (b) All the mass points of the 
disk initially consist of a triangular lattice. When the deformation occurs, it is 
possible that the configuration of mass points (points in figure) locally change 
like this figure. Note that these two configurations are energetically equivalent. 



41 



Fig. 1. authors: Hisao Hayakawa & Hiroto Kuninaka 
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43 




Fig. 3. authors: Hisao Hayakawa & Hiroto Kuninaka 
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